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ABSTRACT 


Grid generation technique using elliptic equations 
is discussed and applied to generate grids on complex 
geometries such as an elliptic cylinder containing a 
circular cylinder and an Aero-Assisted Flight Experiment 
Vehicle. The effect of the control functions on moving 
the grids to the desired locations is also discussed. 

It is hoped to couple this computer code with a flow field 
at some later stage. 


Adaptive grid technique using equidistribution scheme 
is applied to the hypersonic stagnatior point merged layer 
problem which had been analysed earlier on uniform grids by 
Jain and others, A weight function dependent upon a combination 
of the gradient and curvature of the profiles of temperature 
and normal component of velocity is used to place larger 
number of grids near the surface and in the outer shockwave 
like region. Computer code of Jain et al was accordingly 
modified and the resultant code was executed for shuttle 
flight conditions from an altitude of 104.93 Km to 85.74 Km. 

The results on adapted grids are compared with the results 


on uniform grids. It was found that upto a certain stagnation 

Poo 


Reynolds number (Re, 




= 160) corresponding 


to shuttle flight conditions between 92 Km and 85 Km, the 


results with uniform grids and adapted grids are exactly the 
same. The computer took more time to execute the computer code 
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with adapted grids than with uniforTn grids. Within the 
range of stagnation Reynolds number 160 < Re _< 260 
(corresponding to conditions between 92 Km and 85 Km) , 
uniform grids gave converged but inaccurate results while 
adapted grids gave converged and physically plausible 
results. Beyond the value of stagnation Reynolds number 
of 260 and upto 430, computer code with uniform grids 
failed to converge while the computer code with adapted 
grids converged and gave physically plausible results. 

In this manner, adaptive grid technique extended the range 

of applicability of the merged layer approximation, from 
Re = 160 to almost three times its value. Also, 
computations with adapted and uniform grids were carried 
out on a cold wall and an insulated wall conditions with 
varied number of grid points . It was found that the number 
of grid points required with adapted grids to get the same 
accuracy is considerably less than that required with uniform 
grids. Consequently the computer time is considerably 


reduced . 



APPLICATION OF GRID GENERATION TECHNIQUES TO 


(i) Complex Geometric Shapes and 

(ii) Hypersonic Low Density Flow Field in the 
Stagnation Region of a Blunt Body 
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CHAPTER 1 

INTRODUCTION 


I, Grid generation techniques and grid adaptation 
techniques are the modern and powerful tools available 
to solve numerically a variety of problems involving 
complicated body shapes and large solution gradients in flow 
problems . One of the grid generation techniques which 
generates highly acceptable grids is the elliptic equations 
technique. Grid generation using elliptic equations tend to 
produce smooth grids but the grid lines may not be orthogonal 
every where and may be excessively skewed. Orthogonality near 
the surface can be imposed either by specifying the angle 
of the coordinate lines cutting the surface as proposed by 

4 

Thompson or by using additional conditions as proposed 
by Coleman and Haussling^. Skewness can be controlled by 
the proper choice of the control functions in the Poisson 
equations . The advantages of grid generation techniques 
including elliptic differential equation technique are (i) 
the boundaries are represented by coordinate lines and 
(ii) the governing equations are solved in a rectangular 
computational domain irrespective of the shape of the^ ^ ^ ^ 
physical domain. However, additional elliptic partial 
differential equations have to be solved in addition to 
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those governing the fluid flow. 

Winslow^ has initiated the elliptic grid generation 
technique which had been extensively studied by Thompson 
et al^"^. 

In the present investigation, elliptic equations 
(both Laplace and Poisson) are used to generate grids 
about two problems (i) Elliptic cylinder enclosing a 
circular cylinder and (ii) Aero-Assisted Flight Experiment 
Vehicle. For the time being, grids generated here have not 
been coupled with any flow problem. 


II. Due to inadequate resolution of gradients, inaccurate 
solutions are obtained when a flow problem solved on a 
uniform grid, has large solution gradients whose location 
is not known a priori. These gradients can be resolved 
more satisfactorily by making the grids dependent on the 
solution gradients (or /and curvatures) and allowing the 
grid in the physical domain to evolve along with the solution. 
This is called adapting the grid. The grids in the computa- 
tional domain is rectangular, and uniform. The 

governing equations are solved on the uniform grids in the 
computational domain. 

8 '■ 

Brackbill and Saltzman had shown from variational 

principles that the solution error will be a minimum if 

• 2 >■ 

the quantity wJ is equal to the same constant throughout 
the mesh, where w is a weight function which is dependent 
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on the solution gradients (or/and curvatures) and J is 

the Jacobian of the transformation between the physical 

and computational domains. This leads to the equi- 

distribution scheme in the one dimensional case where 

the minimisation of the integral corresponding to the 
2 

quantity wJ is analogous to minimising the energy of 

a system of tension springs, connected in series with. 

each other with spring constants w. 

8 9 10 11 

White , Dwyer et al * , Rai and Anderson , 

12 1314 

Gnoffo , Nakahashi and Deiwert ' and Abolhassani 

15 

et al have developed the equidistribution idea with 

different weight functions and applied to specific 

problems to show the robustness of the method 

in resolving sharp solution gradients . Nakahashi and 

Deiwert^^' ^^liave also developed a procedure for two^^ 

13 

and three dimensional adaptation using torsion spring 

apalogy in addition to tension spring analogy. 

15 

Abolhassani et al have given a form of the weight 
function which is dependent on more than one solution 
variable. 

In hypersonic low density flow, the present 
approach catches the shock wave like structure as part of 
the computational domain. Thus large gradients of the 
flow variables occur near the surface and also in the 
outer region where a shock wave like structure is formed. 
The location of the shock wave in the flow field is not 
known a priori. Hence adaptive grid generation technique 
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seems to be the appropriate technique for solving the 
merged layer problem. 

The following form of the weight function, as 

8 

proposed by Dwyer , is used in this investigation: 


w(C) 


i-M If I 


B 


d^f 


1 2 

dr 


where ^ is the computational coordinate line in the 
computational domain and f is the dependent variable 
which is chosen as the temperature near the surface and 
the normal component of velocity away from the surface. 
Numerical experimentation with various dependent 
variables indicated the above combination along with A = 1 
and 0,05 B ^ 0.4, give the best resolution of 

the flow field. 


Following Nakahashi and Deiwert , the grid interval 
in the physical plane is related to the weight function 
in the computational plane by the following relation. 


Xi = 


w 


(q) 


N 

Z 

i=l 




N 

where L = ^ Ax., N being the total number of grid 

i=l 


points 
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Merged layer problem has been extensively analysed 

and studied on a uniform grid by Jain^^”^^, Jain and 

Adimurthy^^”^^, Kumar and Jain^^'^^, Jain and Kumar^^, 

2 S 2 7 28 

Jain and Saro5 Prabha ' , Jain and Woods and Singh 

29 

and Jain • In all the above investigations, and in 

the current work flow in the merged layer is described 

by using the full Navier-Stokes equations. Levinski 
30 

and Yoshihara had earlier solved the same problem 
using thin shock layer approximation of the NS 
equation. 

With adaptation, it is found that (i) the range of 

applicability -of the merged layer theory can be extended 
to a much denser region of the atmosphere (ii) the 

accuracy of the solution has increased in the range of 
prescribed parameters where uniform grids do not give 
satisfactory results. (iii) Solution with adapted 
grids can be obtained with the same accuracy as the 
solution with uniform grids using lesser number of 
grid points, thereby reducing computer memory and 


time . 
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CHAPTER 2 

GRID GENERATION USING ELLIPTIC 
PARTIAL DIFFERENTIAL EQUATIONS 

2.1 Description 

The first step in solving boundary value problem 
numerically, is to generate suitable grids in the physical 
domain so that the boundary is properly represented. 

Prior to the introduction of boundary fitted coordinate 
system, only Cartesian coordinates were used to solve proolems 
involving arbitrary bodies where the boundary is 

placed between grid points and values of the dependent 
variables and the boundary conditions are interpolated 
onto the boundary. In most boundary value problems 
this interpolation has led to large errors in the solution • 
of the governing equations . Such large errors especially 
occur in the cases of bodies involving large curvatures 
or/and slope discontinuities. The avoidance of interpolation 
is the major advantage of using boundary conforming 
coordinates. 

Transformation relations have to be obtained in 
order to represent the physical domain using boundary 
conforming coordinates (C/Ti) which is already represented 
by the Cartesian coordinates (x,y) . In general, lines 
of constant-^ and constant-n , in the physical domain, are 
non-orthogonal curvilinear coordinate lines. But in the 
computational domain the same lines of constant-? and 
constant-ri are straight lines and the physical domain in 
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the transforined plane (or in the computational plane) 
becomes a rectangle. The intersections of the lines 
of constant- C and constnt-ri give the grid point locations - 
In the computational plane the grid intervals, and 
Af) , are always taken as unity in both the directions , 

C and n . The actual values of the curvilinear 
coordinates ? and ri are immaterial to the subsequent 
use of this coordinate system in the numerical solution 
of partial differential equations, for the grid intervals, 

A? and Ari» simply cancel out of all difference expre- 
ssions for the transformed derivatives . The governing 
equations are always solved in the computational plane 
and the solution in the physical plane is obtained by 
incorporating the metric terms of the transformations in 
the governing equations . 

Some of the restrictions which must be imposed on 
the transformation relations are (i) there should be 
one-to-one correspondence between the physical and 
computational planes, (ii) lines on which one of the 
coordinate is constant, the other must vary monotonically 
(iii) lines of constant-^ (or constant-n ) should not cut 
each other (iv) coordinate lines must be smooth and (v) there 
should not be any void in the physical plane. 

Thus the problem of grid generation boils down to 
finding the transformations between the physical and 
computational domains . The various methods of grid 
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generation differ only by the manner by which these 
transformations are defined. In the case of algebraic 
grid generation the transformation relationship between 
the physical and computational plane is mentioned 
explicitly, either analytically or numerically- The 
initial grids which are used in this investigation to 
generate the solutions of the elliptic equations are 
examples of algebraic grid generation technique where 
the relationship is fixed numerically. 

Consider Fig. (2.1), Fig. (2.1a) shows the physical 
plane where two arbitrary curves and T2 represent 

the inner and outer boundaries. It is required to generate 
a boundary fitted coordinate system in the domain between 
them. For this a cut is made arbitrarily in between 
them at 0=0 (chosen arbitrarily) . The inner body is 
assigned a value of n , say n = 1. The outer body is 
assigned another value of rj , say t)= J. One edge of the 
cut, represented by F^ is assigned a value of say 
5=1 and the other edge of the cut, represented by F^ is 
assigned another value of 5 , say 5 = 1* All the remaining 

values of n ,from ti= 2 to ri = J-1 and 5 torn 5 * .2 to 5* 1-1 
are assigned arbitrarily along the various curves in 
the physical domain confined between F]^ and F 2 as shown 
intFig. 2.1. Thus the physical domain after this transformation 
becomes the rectangular domain in the computational plane. 

■ -i. ■ 
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Conceptually the physical domain can be considered to 
have opened at the cut Q = 0 and then deformed into a 
rectangular domain. This implies that the boundaries 
on which ^ = 1 and ^ = I in the computational plane 
becomes reentrant on each other in the physical plane 
and therefore the boundary conditions need not be and 
should not be specified on them. 

So, it could be seen that boundary conditions are 
specified on all boundaries where it could be specified 
(i,e. on lines with ri = 1 and p = J) . Such a problem 
could obviously be solved using an elliptic system of 
partial differential equations. 

From variational principles, smooth grid lines are 
given by the minimisation of the following integral 

Ig = /[(VC)^ + (Vh)^] dx dy (2.1) 


in two-dimensions . Consider 


I ^ (VC)^ dx dy 
D 


( 2 . 2 ) 


In two-dimensions. 


V? 


i > i ^ 

8x -* 3y 


which yields 
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(Vi) 


2 , ( 3 ^ ^ 2 

^ ^3y ^ 


Hence, equation (2.2) becomes 


~ fi ^ (2.3) 


In general if 


I = // F(x,y,z , dx dy 

D 

then minimisation of I leads to the Euler-Lagrange is equation 


_A' 

3x 


{F } - 


~ — {f T- 
3y ^ 


(2.4) 


where 


3z 

9x 


q = 


M 

3y 


and 


ffpl = V + '■pz li + ^pp If + ^pq |5 

(2.5) 


^{F„} = F„„ + F^., + . Itr + F„„ IB + F, ^ 


ay. q^ 


qy 


qz ax ■ ‘qp ay v ^qq ay 

( 2 . 6 ) 

In equations (2.4), (2.5) and (2.6) suffixes denote 


d i f f er ent i at ion . 
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Applying the Euler-Lagrange equation (2.4) to minimise 
the integral in equation (2.3), one gets 

= 0 (2.7) 

By adopting a similar procedure for the other coordinate 
line h, one gets 

V^n = 0 (2.8) 

Thus, minimisation of the integral (2.1) yields, equations 
(2.7) and (2.8), which are the familiar Laplace equations 
and they are elliptic in nature. It was mentioned earlier 
that the given problem of grid generation in the physical 
plane can be solved by using a system of elliptic equations 
and now it is shown that the particular form of the elliptic 
equations, to obtain smooth grid lines, is the Laplace 
equations or more generally it could be the Poisson 
equations. 

f 

2.2 Mathematical _Developm:eht 

Transforming the two-dimensional, doubly connected 

region D, (physical plane) bomded by two simple, closed 

arbitrary contours, (inner boundary) and T2 

* 

(outer boundary) , onto a rectangular region D (computational 
plane) as illustrated in Pig. (2.1) . is considered here. 




y 




Cb) Tronsformed plane 


FIG. 2-1 FIELD TRANSFORMATION - SINGLE BODY 
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. ★ ★ 

It is required that map onto » ^2 ^2' 

★ ★ 

F- onto F^ and F. onto F.. It is to be noted that 

3344 
★ ★ 

and F^ are required to be constant ri-lines, while 
the arbitrary cut between contours F^^ and F 2 (i.e. 
F^ and F^) becomes oconstant 5-lines. 


If it is required to generate the curvilinear 


coordinates (5 ,ti) in the physical plane, then, the 


elliptic system is of the form 


^xx + ^yy * 

^xx yy Q(5 »ti) 


(2.9) 


with Dirichl dt boundary conditions, one coordinate being 
specified to be equal to a constant on the inner 
boundary and equal to another constant on the outer 
boundary, with the other coordinate varying monotonically 
over the same range around both the inner boundary and 
the outer boundary. Note that in the above equations 
5 and t) are the dependent variables. 

Since the grids in the physical domain is unknown 
and it is desired to perform all numerical computations in 
the uniform rectangular transformed plane, the dependent 
and independent variables must be interchanged in 
equation (2 .9) . This results in the coupled system 
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- 2 6 +• = -J^ [ p(?,ri) + x^ Q(C,ti) ] 

( 2 . 10 ) 

- 26 x^^+ -J^[ y^ p(C,n) + y^ 0(C,n)] 


where 


2 2 2 2 
e = Xj + yj y^, J= Xj y^ _ y^ 


with the transformed boundary conditions 




fl (C, 1 ) >1 

*2 <5- 1 V 


•r? -1 h Tj 


( 2 . 11 ) 


/gi (?. J)1 

( 

\ Q2 ^ # J ^ j 


[E, j ] e r; 


The functions f^, f 2 # g^ and g 2 are specified by the 

known shape of the contours and r 2 and the specified 

distribution of 5 thereon. As noted, boundary data 
are neither required nor allowed along the reentrant 
boundaries r_ and F.. In equation (2.11) J denotes 
the number of grid points in the n -direction. 
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The system given- by the equations (2.10) is a 
quasi-linear elliptic system for the physical coordinate 
functions, x (5,0) and yCC,!!), in the transformed plane. 
Although this system is considerably more complex than 
that given by equation (2.9), the boundary conditions are 
specified on straight boundaries and the coordinate spacing 
in the transformed plane is uniform. The boundary-fitted 
coordinate system generated by the solution to (2.10) has 
a constant n-line coincident with each boundary in the 
physical plane. The ^ = constant lines may be spaced 
as desired around the boundaries, since the assignment 
of the 5-values to the [x,y3 boundary points via the 
functions f^, f 2 , 9^ ^2 arbitrary, (Numerically, 

the discrete boundary values transformed 

to equispaced discrete ~ points on both boundaries) . 

Control of the spacing of the h = constant 

lines and of the incidence angle of the 5* constant 
lines at the boundaries is accomplished by varying the 
functions p(5,n) and Q(5#n) in (2.10). 

2. 3 Control Junctions 

The functions P and Q are called control functions 
as they control the positioning of the curvilinear 
coordinate lines. When both P and Q are zero,, the 
resulting equations are the well known Laplace equations. 
With the Laplace ^uat ions the coordinate lines in the 
physical plane will tend to be equally spaced in the absence 


y 



Q<0 P<0 

FIG. 2-3 


FIG. 2.2 COORDINATE LINES NEAR A (a) CONVEX BOUNDARY 

^ BOUND ARY 

FIG. 2. 3 EFFECT OF P AND Q IN MOVING THE COORDINATE 
^.INES 
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boundary curvature because of the strong smoothing 
effect of the Laplacian. They will be more closely 
spaced over convex boundaries and less so over concave 
boundaries as illustrated in Fig. 2.2. On the convex 
surface. Fig. 2.2a, it could be seen that ri > 0 , 
because of the convex curvature of the lines of constant 
n . Therefore it follows that < 0, and hence 

the spacing between the t] - lines must increase with y. 

The ri -lines will tend to be more closely spaced over 
such a convex boundary segment. For concave segments, 
illustrated in Fig. {2.2b), it could be seen, h „< 0, 
so that riyy > 0, and hence the spacing of the ri-lines 

must decrease outward from this concave boundary. 

Now, suppose that both P and Q are non-zero. In such 
cases, negative values of the function Q will cause the 
n-lines to tend to move in the direction of decreasing p 
while negative values of P will cause C -lines to tend to 
move in the direction of decreasing g . These effects 
are illustrated in Fig. (2.3) for an p-line boundary. 

One particularly effective procedure is to choose P and 
Q as exponential terms and this procedure is followed 
in this work. The coordinates are generated as the 
solution of equation (2.10) with P(C,n) and q(^,ii.) given by 
the following expressions : 



sign (C - exp(-Cj^|^ - C^I) 

b. sign exp(-dj ((C-Cj)^+(r]-Tij)^) 

(2*12a) 

3-^ sign (n-rii) expC-c^ln- ) 

b^. sign(n-nj) ejq) (-dj ( (^-5 j ) ^ + (n-Tijf 

(2.12b) 

where the positive amplitudes and decay factors are not 
necessarily the same in the two equations. Here the 
first terms have the effect of attracting the C ~ constant 
lines to the g * lines in equation (2.12a), and 

attracting i*= constant lines to the ri = lines in 

equation (2.12b) . The second terms cause C * constant 
lines to be attracted to the point (§j^/Ti^) in (2.12a) with 
similar effect on ti= constant lines in (2.12b). 

Moreover, the grids can be made to depend on the 
solution of the governing equations by incorporating the 
solution variables in the functions P and Q. It is 
interesting to note that this becomes some sort of 
adaptive grid technique, in such cases . 


P (?#n) 


s 

i»=l 


m 

Z 

j=l 


Q (§,n) 


n 

- Z 
i*=l 


m 

Z 
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2.4 Nurnerical iTnpleTnentation 

All the derivatives in equation (2.10) are approximated 
by second order central difference expressions 



% 

2 ("kl,j - fk2,j) 

ij 

'll 

i ^^i,Jtl - fi,j-l) 



Sl,j “ ^ ^i,j ^k2,j 



^i,j+l " ^ ^i,j ^i,j-l 


Ck 

4 ^^kl,j + l “ ^kl,j-l ~ ^k2,j + l ^k2,j-l^ 


kl and k2 are described below. 

The resulting set of 21 ( J-1) non-linear difference equations, 
two for each point (i,j) for i = 1,2,3, ....I and j = 2,3,4.... J-1 
are solved by successive accelerated replacement (SAR) method 
which is described in the numerical implementation section 
of the merged layer problem/. to follow. To solve the equations 
by SAR method an initial grid has to be developed in the 
physiGal plane. The generation of the initial grid is described 
in the section on application of the differential equation 
method. 

It is to be noted that both the lines corresponding to 
i = 1 and i = I represent the cut in the physical domain. 

In other words the same line which represents the cut in the 
physical domain is represented by two lines corresponding to 
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i = 1 and i *= I in the computational domain. So when 
calculating the derivatives corresponding to i *= 1 and 
i * I, kl is taken to be 2 and k2 is taken be I-l. For 
all other values of i , kl = i+1 and k2 = i-l. 

Convergence : 

SAR method is applied to one ?-line at a time from 

i s= 1 to i = I . Both the values of x and y are changed 

along a line. When the maximum error among the values 

of X and y along the line is less than the prescribed 
-4 

limit of 10 , then the solution is considered to 

have converged along that line and the next line is taken . 

In this way convergence is attained on all the 5 -lines 

from i =1 to i =1. Now since lines corresponding to i * 1 

and i * I represent the same physical line, the maximum 

error in the values of the difference of x and y for these two 

lines is found. If the error is less than the prescribed limit 
-4 

of 10 than the solution is considered to have converged 
globally. Otherwise the whole cycle is again started with 
i = 1 , until global convergence is attained. If the 
solutions do not converge to the prescribed limit 
in 15 cycles for the elliptical cylinder containing a 
circular cylinder problem and 60 cycles for the AFE vehicle 
problem, but the error decreases continuously, then the 
solution at the end of 15 and 60 cycles are accepted as 
reasonable for the respective cases. 
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2 . 5 Application of the Elliptic Grid Generation Technique to 
Bodies of Complicated Shapes 

(1) Elliptical Cylinder containing a Circular Cylinder: 

This problem in two- dimensions becomes an ellipse 

containing a circle.lt is desired to generate grids ih 

2 „2 

X y 

the domain bounded by an ellipse, —k + ^-9 * 1 , with 

r'^ 

1 2 2 

r^ = 6 units and r 2 = 4 units and a circle (x-a) + 

2 2 

(y_b) = r , with r = 2 units and (a,b) varied among 

( 0 , 0 ), ( 1 , 0 ), ( 0 , 0.5) and ( 1 , 0.5), Figs ^ 2 . 4 ) - ( 2 . 9 ) , by solving 

equations (2.10) with both P and Q as zero and with the 
following boundary conditions: 

(i) On the inner boundary, B =. 1, 

^ r(n= 1) Cos [o(5i)] 

, i = 1,2,3, ......I 

^i,! r(n= 1) Sin[0(5i) ] 

(ii) On the outer boundary, n = J 

^i,J ri(n= J) Cos [©(?,)] 

= , i = .1,2,3, ... . .I 

Sin lo(Ci)] 

where the functions r(T)), .rj^(r)) r 2 (Ti) and ©(g) on the 
boundaries B = 1 and n = J are specified numerically. 
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For example, the boundary values of x and y at (i, 1) 
in the computational plane are (x. and y. ^ respectively 
for all values of i- Similarly, at the point (i, J) they 
are x. and y. ^ respectively for all values of i. The 

X ^ «J 1 ^ u 

physical plane thus gets transformed to a rectangle in 
the computational plane. 


Initial Grid; 


With AO = 15 / 25 grid points are chosen on the 

inner boundary, starting from 0=0° and ending at 
O = 360°, using equations (2.13a). Similarly 25 points 
are chosen on the outer ooundary . 0 is increased in the 

anti -clockwise direction. Each point (x. i -i) on 
the inner boixndary is connected to the corresponding 
point (x. -r » y. t) on the outer boundary by a straight 

X ^ J X ^ u 

line. These lines represent constant-^ lines. Each 

constant-C line is divided into (J-1) equal parts (with 

J = 16) so that there will be J number of grid points 

on each line. The (x,y) coordinate of the j grid point 
til 

on the i grid line is given by 



(J-DXj^j t 


(J-1) 


(j-i)yi,j + 
(J-1) 


, j = 2,3,4 (J-1) 

(2.14a) 

, j = 2,3,4 (J-1) 


(2.14b) 


Now by keeping J = constant and varying i, all the 
points (xj_^j ' ^i,j^ are joined with each other. For each 
value of j from j = 2 to j = J-l, this is repeated. These 
lines are the constant-rj lines. 









FIGURE £.8 Centre of the Circle at C 1 ^ 0 . 5 ) “In it La 1 Gr 
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Results : 

Initial grids are shown in Fig. (2.4) and Fig. (2.8) 
for the cases when the centre of the circle, (a,b) lies 
at (0.0) and (1, 0.5). Final grids are shown in Figs. (2.5), 
(2.6), (2.7) and (2.9), when the centre (a,b) lies at 
(0.0), (1.0), (0,0.5) and (1,0.5) respectively. By comparing 

the final grids with the initial grids it could be seen 
that the n-lines near the inner boundary are almost circular 
in shape and those near the outer boundary are nearly ellip- 
tical in shape. This is due to the inherent smoothing nature 
of the Laplace equations . CPU-time taken in the ND-500 system 
of the CAD-center is about a minute. 

' f 

( 2 ) Aero-Assisted Flight Experiment Vehicle ; 

In order to investigate the flow characteristics 
around a reentering vehicle, NASA is planning to carry out 
exhaustive flight experiments using a vehicle similar to 
the one shown in Fig. 2.10. This vehicle essentially 
consists of a raked cone followed by a circular cylinder. • 

It is a fairjy complicated shape in that it is a combination 
of several bodies. It has sharp corners and a concave region. 
The inner and outer boundaries are generated as follows. 

Inner Boundary : 

Inner boundary is constructed with respect to the axes 
ox and oy in Fig. 2.10. BCD is an arc of an ellipse 
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X V T 

— 2 + = 1 with = 1.0 units and r 2 = 0.6 

^ 1 . .^2 , , 
units with centre at O. The coordinates of 

the point A, where the raked surface cuts the x-axis is 

given by (-0,6, 0,0). The equation of a straight line 

which passes through the point A and makes an angle Q with the 

x-axis is given by y = m(0) (x - x^) where m(0) * tanO. 

The intersection of this line with the prescribed 

value of O ( ©= Oq = 70° or 90°) , with the ellipse gives 

the points B and D. The x value of the points B and 

D is given by 


Tn(0) 


X. 


X = 


1^2 + rf lm(0) 1 


2 - 




m 


(0)1^. (r| 






^2 + ri {m(0) 


with the (+) sign for B and (-) sign for D. These values of x 
are substituted in the equation of the straight line to 

get the corresponding value of y. 



The equation of the line HG is y = 0.2 x AA and the inter- 
section of this line with Sb gives the point H. 

■ • ■ ■ 

■ •' ' 0 • 2 X AA. 
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Another line EF, y = -0.2 x AA is also drawn and 
the intersection of this line with ^ gives the point E. 

Yg = -0 .2 X AA 

^A 

The intersections of the line x = 0 with the lines 
y “ Yjj and y * y^ gives the coordinates of the point G 
and F respectively. Thus, x^ = 0.0, y^ = y^^, Xp = 0.0 

and Yp = Yp. 

It is decided to distribute N1 (= 4), N2 (=8), N3 = (5) 

N4(=9), N5(=5), N6 (=8) and N7(=3) grid points on the segments 

5g, GH, Sb, BCD, DE, EF and FO respectively. The following 

convention is introduced here to make the discussion simpler 

and easier to understand. N12 = N1 + N2 , N13 = N1+N2+N3 

and SO on and-Nl? = N1+N2+. . .+N7. Thus point o is the 

first grid point, point G is the N1 grid point, point H is 

the N12^^ grid point, and so on, and point F is the N16^^ 

til 

grid point and again point O is the N17 grid point. 

The quantities ^Yq* and Axp are defined as follows: 

LYq « OG/(Nl - 1) 

AXjj = GH/N2 
AXp = EF/N6 

They will be used while defining the boundary conditions. 
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Outer Boundary : 

Outer boundary is constructed with respect to the XO "y 

axes- 

The outer boundary consists of a parabola IKM and 
a straight line TH. The equation of the parabola is given 
by 

2 

y = ax + c (2.15) 

The values of the constant a and c are found by 
specifying two points K and J on the parabola. The 
coordinate Y^) of the point K is (0.5, 0) which gives 

from equation (2.15) c =-aXj^. Coordinate (Xj, Yj) of 
the point J is prescribed as follows: 


= ’'b 

BJ is perpendicular to x-axis and is assigned a value 0.6. 
Substituting the values of Xj, Yj and C in equation (2.15) 
the value of a is found to be 

a = Yj/(Xj - X^). 

To distribute the grid points, the outer boundary 
is divided into segments 0*1, IJ*, J KL, LM and MO . The 
points J* and L are obtained by the intersection of a^ ^ 
straight line, making an angle 80° with the x-axis and passing 
through the point A, with the parabola. 
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The x-coordinates of the points J' and L is 
given by 


X = (- 


+ 


2 lm(0) 





4 m(G) 


(c+a X^) 


2 lni(o)l 


+ sign is taken for the point J' and -sign for the point L. 
Corresponding Y values are found from the equation 
Y = m(0) (X - X^) . 

The number of points that are distributed on the 

segments 5^, J*lE, LM and MO* are (= 8 ), M 2 (* 9), 

M^ (= N4 = 9), M 4 (=9) and M^ (=7) respectively. Thus, 

th 

point O* is the first grid point, I is the Ml grid 

tltl til 

point, J* is the M12 grid point, L is the M13 grid 

th 

point M is the M14 grid point and again O' is the 
til 

M15 grid point. The quantities Ay^, AXj,, AXj^ 
are defined as follows and are used while defining the 
boundary conditions. 


AYj = 

Yj/(M1 

- 1) 

AXj , » 

(X^. - 

Xj)/M2 

AXl = 

(XqI - 

x^)MA 
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Boundarv._Conditions : 

The distance KC in Fig. (2.10) is defined to be 0.5. 
This makes the transformation relationship between (x^y)-axes 
and {X,Y)-axes to be 

X = X + o"o = X + 2.0 
Y = y . 

The boundary conditions on the inner and outer boundaries 
are defined with respect to (X,Y)-axes and the subsequent 
generation of the grids is also with respect to (X,Y)-axes 
only. 

(i) Inner boundary i.e. ri = 1 

(a) Segment OG. 

X. . =2.0 with = i = 1,2, ....Nl 

X ^ X X 

Y^ = (I-DAYq with = i = 1,2. ...Nl 

(b) Segment GH. 



= (i-Nl) AXjj 

with 

= i = Nl+l,.. .,N12 


= ^N1 

with 

K i ^N1 +1, ...,N12 


(c) Segment HB 

(i^Nl2)Xg + (N13.-i)Xjj 
Xi,i - ' - 5l 

= ,m2+l, . . .,N13 
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Y 


irl 


(i-N12) Yg + (N13-i)Yjj 
___ 


with 


i *N12+l5--N13 


(d) Segment BCD 


©j^ = Qq + {i-N13 ) AO ' where AQ = 20 




rl +fm(e^)p X 


rir 2 jr 2 +|m(Q ) P (rf-x|) 


2 2 
^2 ^ 


(Ol)| 


2 2 ! /„x r 2 


^ 2+^1 im 


(O) 1 


with * i « N13+1, ,N14 


-ve sign is taken if Tr/2 — sign 

is taken otherwise. 

Y. = M(0.) (X. - - X ) with = i « N13+1 N14 

X ^ X X X ^ X A, 


(e) Segment DE. 

(i - N13 )Xe + (N14-i)Xg 


==1.1 


Y. , 

1,1 


(N5-1) 

(i - Nl3)Yg + (N14-i)Yjj 


(N5-1) 


with ^~i = N14+1,-..N15 


with 5i=i = N14+1,...N15 


(f) Segment EF. 

X . - = (i-Nl5) AXp with * i * N15+1, . . .,N16 


^i,l “ "^NIS 


with * i * N15+1, . . .N16 
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(g) Segment FO. 


""i,! = "^Nie q = i = N16 + 1...N17 

= (i-N16)AYQ with = i = N16 + 1 ... N17 


(ii) Outer Boundary i.e. = J. 

I 

(a) Segment 01. 

X. _ = X , = 4.0 with C,- = i = 1 Ml 

1 / J o' 

Y. . = (i-DAY:^ with = i = 1...M1 

1 / J ± 1 


(b) Segment IJ 


X. 


. _ = (I-MDAXt, with = Ml+1, ...,M12-1 

1 / u ^ J- 


Y. 


. _ = / a X. T + <^ with C- = Ml+1, . . . ,M12-1 . 

1 , J 1 / >J ^ 


(c) Segment J KL. 


= 80° + (i-Ml2)AO where A© = 20° 


. X. 




a^ + 4 1 m (©. ) I ^ (c+a X, ) 
( . ^ — .2 + X ^) ± — 


•2|m (0^) j 


2lm{©.) 1 


with - i = M12, ...M13 
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{-) sign is taken if 'n/l < < 3tt/2 and (+) sign is 

taken otherwise, 

^i,J = with = i = m12 M13 

(d) Segment LM. 

^i,J "" (i-Ml3)AXj^ with = i = M13+1,...M14 
= Va j + c with = i + M13+1,...M14 

I 

(e) Segment MO . 

^i,J " ^i,M14 Ci = i = M14+1 M15 

= (i-Ml4)AYj with = i = M14+1, M15. 


Initial Grid and Solution Procedure : 

Points (xjj^ j^) on the inner boundary are 

joined to points j/ Y^ j) on the outer boundary, with 
J =36, by straight lines- These are the constant-^ lines. 
Each constant-^ line is divided into ( lP- 1) equal parts and 
the constant-n lines are constructed as described in the last 
example. Using this grid, SAR method is used to solve the 
Poisson equations (2. 10) with P= 0 and the following 
cases of Q: \ ' V' 
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(i) Q = 0, This gives Laplace equations. 

(ii) Q as given by equation (2.12b) with b^ = 10.0 
and = 0.025. Constant-ri lines are to be attracted to 
the point H in Fig. (2.10). 

Discussion of the 'Results : 

The initial grid, solution of Laplace equations and 
solution of Poisson equations around the AFE vehicle are 
shown in Figs. (2.11), (2.12), and (2.13) respectively. 

In these figures and in Fig. (2.10) the outer boundary IKM 
(Fig. 2.10) is a parabola. But there seems to be a sharp corner 
at point K. . This is due to the fact that the adjacent 
grid points are connected by straight lines, instead 
of smooth lines, by the computer. 

It could be seen that a marked difference exists 
between the initial grid and the other two. The grid 
lines are smooth in both the cases of Laplace solution 
and Poisson solution. In both the Figs. (2.12) and 
(2.13) grid points are concentrated near the surface 
where the gradients of the flow variables are supposed to 
be large and sparsely distributed away from the 
surface. ■ 

Comparing Figs. (2.12) and (2.13) it could be seen 
that the first constant-Ti line away from the inner body 
is almost the same in both the cases but other constant -q 







FIGURE 2.12 ftero-Assisted Flight Experiment UehicLe 







FIGURE 2,'iJ ftero-ftss isted Flight Experiment 
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lines are attracted towards the point H, as seen in 
Fig. (2.13). The inability of the first constant -H 
line away from the inner body to move closer to the 
inner body, in the region near the point H, may be 
due to the sharp concave corner at point H of the inner 
body or due to the fact that the choice of the control 
function Q in Poisson equation (2.12b) needs further 
modification by numerical experiments. More constant-n 
lines can be attracted towards the point H, either 
by taking a smoothly varying concave surface around 
the point H, instead of the present sharp corner, or 
by making the control functions P and Q in the Poisson ■ 
equations (2.11) dependent on the flow variables which 
are obtained as the solutions of the governing equations. 
For similar reasons, it is observed that there are few grid 
points near the concave corners I and M of the outer 
boundary . 
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CHAPTER 3 

ADAPTATIVE GRID GENERATION TECHNIQUE 

3.1 Description 

In Chapter 2/ the problem of errors arising due to 
inaccurate representation of the boundaries of bodies 
of arbitrary shapes was addressed to, by developing a 
body'^fitted coordinate system. In the present chapter, 
the problem of errors, in the numerical solution of 
flow fields, arising due to inaccurate resolution of 
the gradients and/or curvatures of the flow profiles is 
addressed to. Generally, such large gradients and/or 
curvatures occur near the surface and in certain cases 
they in addition occur somewhere in the flow field 
whose location is not known a priori. These gradients 
and/or curvatures of the flow profiles can be resolved 
properly by moving more number of grid points to regions 
where gradients and/or curvatures are large from regions 
where they are small. In an iterative procedure of 
solving the flow problem, the grids are moved to desired 
location as the solution develops from some initial 
stato. 

As mentioned in the last chapter one way of moving 
the grid points is to make the control functions P and Q 
dependent on the gradients and/or , curvatures . But this 
approach is not quite attractive because it is an indirect 
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way of moving the grid points and it is very difficult 
to choose the particular form of the control functions which 
will move the points to the desired locations. 

A more direct method would be to relate the gradients 
and/or curvatures of the solutions of the flow problem 
to the grid intervals in the physical domain in such a way 
that the grid intervals become small at locations of large 
gradients and/or curvatures. Such a method of moving the 
grid points is called the adaptive grid technique. 

Since the grid points in the physical domain are 
moved, the grid intervals become uneven and the coordinate 
lines may become non-orthogonal . For these reasons the 
physical domain is again transformed to the computational 
rectangular domain as in the last chapter. The initial 
transformation is fixed numerically in most of the cases. 

In the present investigation adaptive grid technique 
is applied to a one-dimensional problem. The physical 
coordinate is represented by n and the computational 
coordinate is by C » Both n and C vary from zero to 
unity with the grid size in the computational domain being 
A5= 1/NDIV where NDIV is the total number of grid intervals. The 
grids in the physical plane are moved depending on gradients 
and curvatures or gradients alone. 
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3.2 Mathematical Development 


7 


As shewn by Brockbill and Saltzmann , for one- 
dimensional adaptation*, the total error is reduced when 
the grid points are distributed such that 


X 

^i+1 

w (x) dx = constant (3.1) 

^i 

for all i* where w(x) is a positive weight function. 
Usually the gradients of the evolving solution are taken 
as the weight function as they are related to the error 
in the solution* 


To maintain a uniform grid in the C-space equation 
(3.1) implies that 

X w(5) = constant (3.2) 


where x^ is the metric coefficient and corresponds to 
the ratio of arc lengths in the physical and computational 
space. But equation (3.2) is also the Euler-Lagrange 
equation for minimising the integral 


^w 


I 



(3.3) 


physically the minimisation of the above integral represents 

the minimisation of the energy of a system of tension 

springs with spring constants w(5) between each pair of grid 

points. The constant in equation (3.2) is evaluated 

14 

according to Nakahashi and Deiwert as follows . 



35 


Let the constant be denoted by c. Rewriting equation 


(3.2) , 


dx. = 


c d g 
w(g) 


(3.4) 


Integrating over all the grid intervals. 


/ <^^1 


1 

f 


dg 

w(g) 


1 .e . 


1 

/ 

o 


_1L- 

w(g) 


(3.5) 


where L is the total length in the physical plane, 
Rewriting equation (3.5), . 


c » L/ / 
o 


IL 


w(g) 


(3.6) 


Also, if there are N grid intervals then Ndg = 1, 


which implies dg 


N 


Putting (3.6) and (3.7) in (3.4), 


or 


dx^ * l/[n w(gi) / ^) ] 

= !■/[ w<Ei) E 1 

1 = 1 ■■ ^ ^ 


(3.7) 


(3.8) 


Equation (3.8) Is the equation used in this work for 
grid adaptation. The form of the weight function is taken 


as 
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w(Cj_) * 1 + A I If I . + B 1^1 (3.9) 

i dC i 

with A to be one in all cases and B varying from 0 to 
0.4. f is chosen to be a combination of temperature 
and normal .... component of velocity. The reasons 

for choosing this particular combination of the dependent 
variable in the weight function are stated in the section 
entitled "choice of the dependent variable in the weight 
function w" in Chapter 4. 



37 


CHAPTER 4 

HYPERSONIC MERGE_D LAYER FLOW NEAR THE 
STAGNATION POINT OF A BLUNT BODY 

4.1 Introduction 

When a space vehicle reenters the atmosphere at about 
120 km altitude, its velocity is of the order of 6 to 10 km/s. 

— 9 

The ambient conditions are (i) the density of the order 10 

3 

kg/m and the molecular mean free path of the order of Im. 

The boundary layer thickness is inversely proportional 
to the square -root of Reynolds number and the shock thickness 
is directly proportional to the freestream molecular mean 
free path. The Reynolds number is very small since the density 
is very small eventhough the flight speed is large. This means 
that the thickness of the boundary layer is very large and 
therefore boundary layer approximations are not valid. The 
shock wave is very thick since the mean free path is large and 
therefore Rankine-Hugoniot relations are not valid. Moreover 
owing to their large thicknesses both the viscous region near 
the surface and the shock wave like region merge into each other 
and this makes the entire disturbed zone viscous. The whole 
disturbed zone is called the merged layer. 

Particles moving in the merged layer do not suffer enough 
collisions with other particles and by the time they reach 
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ths body / tlisy rotain sttough mo^nantuni to slida ovsr the body 

with a finite velocity. This velocity is called the slip 

velocity. Similarly the temperature at the base of the 

viscous layer is different from that of the body. This is 

called the temperature jump. Using kinetic theory approach, 

the expression for velocity and temperature jump are given 
31 

by Kennard as follows: 


% “ (|f)b 


“ ^b ^ay^b 


where u and T are the tangential component of velocity and 
temperature respectively and y is the coordinate normal to the 
body. The subscripts b and w denote the base of the viscous 
layer and the wall respectively. X is the mean free path. 

The flow in the merged layer is described by the full 
Navier-Stokes equations. Local similar solutions valid in 
the stagnation zone of a spherical surface have been 
obtained. The resulting equations are integrated from the 
surface #with slip and temperature jump conditions / to the 
free stream so as to catch the shock wave like structure 
in the outer region as part of the computational domain. 

Merged layer problem has been extensively analysed and 

" 1 Q 20—22 

studied by Jain , Jain and Adimurthi , Kumar and 

Jain^^'^^, Jain and Kumar , Jain and Saroj Prabha'^ ' , Jain 

' ' ^ ' 28 ' ' ' ' ' ' 29 

and Woods and Singh and Jain , using full Navier-Stokes 



equations with slip and temperature jump surface conditions^ 
on uniform grids. 


In the present investigations, computations have 
been carried out with uniform and adapted grids for shuttle 
flight conditions. It had been found that with adaptation 
(i) the accuracy of the results increased, (ii) the range of app- 
licaoilxty of the merged layer approach had been extended 
to a higher range of Reynolds number and (iii) the computer 
time was saved. 


4.2 Mathematical Development 

The full Navier-Stokes equations in spherical polar 
coordinates, with axial symmetry are the following: 


Continuity equation : 


(pv) 


2Pv 


. Pu Pu Goto 

-j- — — ■■ 

r r 


= 0 


(4.1) 


Radial momentum equation ; 


Pr + P (V S V - I-) 


u 


= 2[uv^ - ^ 2v + u^ + u Cot O)]^ 

. l / rr ^ (2v - 

+ r rn >0 r r [ 


+ J {Cot© [r (|)j. + ^ ] - p (V + u Goto) } 

(4.2) 


Tra nsverse moment uTn ecfuation: 


O 


+ P(v u + ^ u. + 

r r Q r ' 


- 2 . r ^ j. \ ■ , 

~ r L r ^ 0 - 3]F (r + 2 v + Uq + u CotO) 

\ 

+ {w |r + -e, ^ 


+ ^ CotO (u - u Cot©) 

r ^ 


(4.3) 


Energy equation ; 


P (v + p h^) 


= V Pr + J Pq + ^ ^ri^hr) r+ r +^^0 


+ "^ (Uq + v)^ + ■^ (v + u CotO)^+(r (p)j. + 


2 P 


3 r' 


[ r + Vq + 2v + u Goto] 


(4.4) 


Equation of state ; 


p = p RT 


(4.5) 


Visoosi tv-temperature relation i 
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Here, r and Q represent the radial distance and vectorial 
angle respectively (See Fig. 4.1); u and v are the velocity 
coniponents in increasing r and Q directions respectively; p , p 
and h are the density# pressure and specific enthalpy respect- 
ively ; P is the coefficient of viscosity; R is the gas constant 
Pr is the Prandtl number; C is the specific heat at constant 
pressure and Y is the ratio of specific heats. Suffixes r 
and 0- denote the derivatives of the dependent variables with 
respect to *r* and '0* respectively. 


Boundary conditions are: 

At the outer edge of the merged layer, r » 


u = SinO 

00 

v = -Vjjg CosO 

p = p (4.7) 

00 

p = p<» 

T = T 

oo 

At the surface of the body, r = rg, slip velocity and 
temperature jump boundary conditions are 





^w ^ a ' W+1 ' 2 Lpr Pa ^ar'Jr 


B 


(4.8a, b) 

(4.9) 


where o and e are the momentum and thermal accommodation 

coefficients and T„ is the prescribed wall temperature. 

. ' w 



Outer shock wave like 



p,P,T = Density^ Pressure and Temperature respectively 

v,u = Velocity components in r ,© 

direction respectively 
^ s Viscosity 

Subscript *oo’ denotes freestream conditions 


FIG. 4.1 




42 


For the adiabatic wall case. 


^dr’r^ = 0 


( 4 . 10 ) 


For the noslip conditions, at r = r 


B, 


u = V = 0 ; T = or (|J) = 0 


Non-diniensionalization Proced'ure 

The non-dimensional variables are defined as follows: 



V 

V 


P 


P V 

^00 CO 



y 

y 

Ooo 




Wo.) 


where the suffix denotes the value of the variable at 

the freestream location and is the freestream stagnation 

temperature. 

The governing equations (4.1) to (4.6) and the 
boundary conditions (4.7) to (4.10) in the non— dimensional 
form become 

(-5, ^ (|u) 9 ^ £u_soto _ „ 

f f XT ■ , " ' f 



p_ + p (vv_ + 7 Vq 

r r x 


u 


Re, 


r y V - (r V + 2v + Uq + u CotO)] _ 

r 3r r ^ r 


O 


r r 


Re^r r 
o 




—{Cot© [r (J) + -^ - - (v + u Got©)} 

L r - - - 


r r 


(4.12) 


u - 

r r 


'© - / — u - , uv. 

— + p (vu + — u_ + — ) 

- w 


2 r li /- j. _ X (r- V + 2v + + u Cot©}J 


- [ ^ (S^ + 9) - - 


Re^r r 


3r 


+ rJ- ^ 


r r 


le,, ^ 3j^r i (a)_ + ^ 

iE ^ 3c R© IT ^ ^ 

o 


3^1 


+ Cot© (Uq - u Cot©) 


(4.13) 


Re^r 



p [A ( V T_ + = T )] . * V p + — 3 

r r f r 


+ 

1 1 A 

Pr Re^ j 

{ (ir T_ 

y Tq 

) _ + + y 

T 

= . y 0 

T + _ 

Goto } 


r 

r r 

r r 


+ 

5 [2tf ^ 
r 

+ -2 
r 

(Uq + v) + -32 

- - 2 
(v+uCotQ) 

— . V ^ 

r=— /ii% . ©T ^ 

+Ir(f) + ri 

r r 


2 u r = 

v_ + V 
r . 

+ 2v + u Goto] ^ 


(4.14) 


-2 ^ ^ 

3Re„ r 


p 

= X-Zi J[ p 

y 

T 



(4.15) 

j: 


= Jt~ 




(4.16) 

where 

A 

P O 00 

- 1 4. 

1 


(4.17) 

A 

2 ■ 
V 

00 

" 2 ■*■ 

(Y- DM^ 



K['©ir0^ ■ '^jp^ 

Yand Pr 

are treated as 

constants. 


Boundary conditions are: 





At f = , r^ 






5 = 

SinO 






- CosO 





p = 

1; f = 

1 - ^ ^ ' 

y^l 

(4.18) 
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p-[A ( V T_ + - f )] . V p_ + 

r r r 


u p 


0 


+ 


1 A 


Pr Re 




o r 


T, 

{ { ;ir T_) _ + + P T + ^ ^ Goto } 


r r 


+ ^[2vf + ZI- (Uq + v)-^ + -^ (v+uCotO)^ +[?{i-)_ + 


r.A2 . 2 


r r 


^ ^_2 I r v_ + -^ + 2v + u Cotol 


SRe^ r r 


(4.14) 


_ Y — 1 — 

P = A p T 


(4.15) 


(4.16) 


where 


A 


P o 1. . 

2 2 
V 
00 


1 


(Y- 1)mJ 


(4.17) 


Here, C_ , Y and Pr are treated as constants , 


Boundary conditions are; 
At'/f = r« 


u = SinO 
V = - CosO 
p = 1; T “ ;1 


'^.•1 

2 A 


P = 


YM„ 


(4.18) 



Slip velocity is 




TTY 


" V 'pyr 


(- 


y 8u 


_)_ 
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(4.19) 


and 

V = 0 

Temperature jump conditions is 


T = T + ) 

^ -^w ' a ^ W+1^ V 2 




Pr Re 
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where 


(4.20) 


(-^)i 

3r r = 1 
(4.21) 
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For an adiabatic wall case. 


{-?) 

3r r=l 


= 0 


(4.22) 


At r =1 and for noslip case 


V = u = 0 and f = T^ or 


3f 

3f 


= 0 


(4.23) 


The r-coordinate ' in the physical plane is normalised 
to the 0 -coordinate which is also in the physical plane by 
the following transformations; 
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where is the total merged layer thickness and n varies 

from zero to unity. This gives 


^ = 1- ^ 1 9 

9r 9ri * 9r n^. * 3ri 

eo 


( 4 . 25 ) 


Method of solution 

The flow and thermodynamic variables are eijpanded 
around the axis of symmetry as follows: 


u (n, 0) = 

Uo(n) 

Sin© 

(4.26a) 

v(n, 0) = 

0 

l> 

Cos© 

(4.26b) 

p(n,©) 

o 

ICL 


(4.26c) 

f(ri, 0) = 

T„(n) 


(4.26d) 

P (n,o) * 

Po'”’ 

- ... 

+ P 2 Sin © 

(4.26e) 

y (n,©) 

5.(n) 


(4.26f) 


The form of various expansions is suggested by the 

, ' ■ , j 

boundary conditions at the outer edge of the merged layer, 

symmetry of the flow around the stagnation line and the fact that 

. ' 2 ' 

pressure inside the merged layer varies as (Sin O) • The later 
statement is conjectured from the following considerations 
(i) Newtonian formula suggests that pressure on the surface 
of a blunt body varies as . Sin^e. Also (ii) pressure behind 
the shock when it is fbrroed varies as Sin^O. This condition 
necessitates including the higher order term p2(Ti)Sin O as 
part of the leading term approximations. For this reason, 
leading term form an indented series. Such considerations are 
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not available to find correct expansion for density. Within 
Newtonian approxi’nation, density is infinite in the thin 
lay©r near the body while Rankine— Hugoniot relation gives 
finite value of density behind the shock when the Mach nu'niber 
is large . These considerations do not lead us to suggest a 
correct expansion for density variation. As such, in the 
present analysis, density variation near the surface bears 
some error. 

The coefficients in the series expansions and their 
derivatives are assumed to be of the order of unity. These 
expansions (4.26a) and (4.26f) are substituted in governing 
equations (4.11) to (4.16). Collecting like powers of O 
the following set of differential equations are obtained: 


1 + n 


eo 


(o V ) + 2p (u +. V ) 

n _ '^o o' ^o o o' 


(4.27) 
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(4.28) 
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(4.32) 
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In the above equation s<a prime over a variable denote 
the derivative of the variable with respect to r\ . These 
equations are strictly valid at the stagnation line only. 


The boundary conditions become. 


at the outer edge of the merged layer, = 1, 


u. 


1 7Vo = 


-1 ; = 1? T 

o o 


1 - 


1 

2A 


= 0 (4.33) 


at the surface, ri= 0 


u. 


/IFT" 2-0 ^ 

2 a Re 


T 

o Ooo 
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eo 


p T 
•^o o 


Q\ 

dn n=o 

(4.34) 


(4.35) 


T = f + (^) ^ Sh" 


a ' Y+1 Pr Re^ 


n 


eo 


dT. 




- =1/2 dn 'ti= 0 

o 


(4 .36) 


For an adiabatic wall case 


.dT. 

^dR' n=o 


= 0 


(4.37) 


For noslip condition , 


"o = ° 


dT, 


To = T^ or ^ 


= 0 at h = 


(4.38) 


0 


(4.39) 
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Note that both the r-line of the ri-line lie in the 
physical plane and the linear relationship between them is 
given by equation 4.24. When the grids are adapted, grid 
points in the physical plane become unevenly spaced. But 
grid points are uniformly spaced in the computational plane, 
by construction. Hence, the governing equations are solved 
in the computational plane. 


Let 5 be the computational line. The relationship 
between f and 5 or 0 and ? are not known, a priori. 

Let, 


U n) 


(4.40) 
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dn 
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(4.41) 
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^ (^) = 
dn 'dn^ 
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r\ 




n 


n 


££. ^ 

3 d^ 


n^ d? 


(4.42) 


Substituting equations (4.41) and (4.42) in equations 
(4.27) to (4.32), the following equations in the computationai 
plane are obtained: 
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(4.43) 
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V, 


T 


(4.49) 


Note that in the equations (4.43) to (4.49) dashes denote 
derivativi^ with respect to ? Shd all the variables, v^» P^* 

Tq» Pq» P 2 » iIq and n are now functions of C . Ihe ccrtputational 
coordinate 5 is also varied froni zero to unity. 
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Boundary conditions 


At C * 1, 


% * Vq = 


-1; p »= 1; T = 1 
o o 


; Po = 0 


(4.50) 


At ? = 0, 
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Vo = 0 


(4.52) 
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For adiabatic wall case 


?= 0 


(4.54) 


For noslip case, at ? * , 0 


io = Vo = 0 


(4.55) 


T^= T,, or 
o w 


(4.56) 



4 , 3 Nunnerical Methot^ 


The governing equations consist of three second 
order differential equations narnely equations (4.44), 
(4.45) and (4.46) for the variables v^, u^ and respect- 
ively . These equations niay be written in the forni 


" _ • _ • 
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%# ) * 0 

(4.57) 

^o 
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_ 1 

* ••••} * 0 

(4.58) 

To 

II 



^o^ K* * 0 

(4.59) 


and the rest are first order equations. 

It is desired to obtain the solution in the range 
Co# d in the computational domains, which is divided 
into equal parts, each of length, A = 


NDIV 


. Thus, the 


total number of grid points is NDIV+1. In equations (4.57) 
to (4.59), the first and second order derivatives are replaced 
by 


N 


^N+1 ^N-1 

2 a 


(4.60) 


^N+1 " ^N-1 




(4.61) 


where Xjj represents the value of any one of the variables 


,th , 


Uo# v^, Tq, Pq, P 2 * Pq or ri at the N 'mesh point. 

Here, N *= 1 represents the surface and N = NDIV + 1 represents 
the free stream location. That is, N = 1 corresponds 'to 
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5 = 0 and N = NDIV + I' corresponds to C = 1 . The difference 
equations so obtained for each grid point are solved by 
SAR Tiethod . The first order equations (4.43) and (4.48) 
are solved by trapezoidal, Simpson and Weddle's rules as 
in section 4.4. 

Method of Successive Accelerated Replacement : 

Equations (4.57) to (4.59) form three sets of algebraic 
equations each set containing (NDIV - 1) number of equations 
(N varies from 2 to NDIV) . The boundary conditions for 
these three variables for N = 1 and N = NDIV + 1 are as 
follows : 


For the slip case 
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(4.63) 
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(4.65) 
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For noslip case. 



U. 


NDIV+1 


1 * V, 


NDIV+1 


■If T, 


NDIV+1 


(4.66) 


1 _ 

2A 


(4.67) 


Each set of (NDIV-1) number of equations is solved by 

the method of Successive Accelerated Replacement (SAR) . 

32 

This method was originally developed by Lieberstein for 

the iterative solution of nonlinear algebraic equations. 

33 

SAR method is a generalization of Young's method of successive 

overrelation of the iterative solution of linear algebraic 
-34 ’ 

equations. Lew first applied this method to solve the 

Blasius equation and obtained 1.7 as the optimum value of 

35 

the acceleration parameter. Later, Dellinger 

applied this method to non-equilibrium viscous layer stagna- 

20 

tion point flows. Jain and Ad imur thy applied this method 

to hypersonic merged layer flow in the stagnation region 
and found that this method gave reasonably accurate results. 

In this method, corrections applied to the values of 
the variables at each mesh point are controlled by the 
acceleration factors which prevent the iterative scheme from 
diverging. This scheme is successive replacement in the sense 
that new variables are made use of as soon as they are available 
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This method decouples the system of algebraic equations. In 
each iteration, it involves the use of one particular equation 
only for the correction of one unknown variable in that 
equation. Also in each iteration the boundary conditions 
are always satisfied. This method is economic in computer 
time as well. Due to all these points the method of Successive 
Accelerated Replacement is applied to the system of equations 
(4.57)tto (4.59). 


Now, the correction procedure for one of the variables, 
say Vq is described here. The other two variables u^ and T^ 

governed by the second order equations will be corrected in 

‘tin "tin 

the same way. The (K+1) approximation to v^ at the N 

mesh point is calculated from 
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The number W is called the acceleration factor, 
^o 


It is calculated as follows : 

First it is required that 
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where e is a prescribed' small parameter- 
Let, 




e 





( 4 . 70 ) 


If o < 1 then W is taken to be and if a„ > l 

o o ”o 

then is taken to be unity. 

This procedure ensures that the correction applied 
at each iteration is never greater than e times the value 
of the variable at that point in the previous iteration. 
The correction starts from the point N = 2 and continues 
upto N = NDIV. 


4 .4 Solution Procedure 

To start with, an initial profile for the flow 
variables satisfying the boundary conditions is chosen. 

The value of is prescribed initially and then corrected 

in successive iteration. Following equation (4.69), the 
corrections are first made to the variables v^, u^ and 
Tg in equations (4.57) to (4.59), in that order, for grid 
points N = 2 to N = NDIV. Then p^ and are calculated 

from the equation of state, equation (4.47) and the 
viscosity-temperature relation, equation (4.49) respectively - 
pQ and P 2 are obtained from the continuity = equation, 
equation (4.43) and the equation (4.48) respectively. Since 



the values of and P 2 are known at the free stream 

end ( Pq =1.0 and = 0.0) but unknown 

NDIV+1 '^NDIV+1 

at the body surface » the integration starts from the 

free stream. 


Now, the procedure for calculating and P 2 will 

be described. At the last but one point from the free 

stream, the values of and P 2 are found by Trapezoidal 

rule. The values at the next four points are calculated 

til 

by using the Simpson's 3/8 rule. At the remaining 
points except the point on the boundary surface, Weddle's 
rule is used to evaluate P 2 and P^ 
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Weddle's rule is used to find the value of P 2 at the surface 
also. 

Due to a singularity in equation (4.71) at the surface 

(• * V =s 0), the density at the surface has a finite 
'°1 

value only if 

V* = -2r\r nec^o ' (4-.73) 

for no-slip case. 


0 
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This condition is used to improve the value of 
second grid point. 


v^ at the 


If equation (4.71) is to be used to evaluate the density 
at the surface, L'Hospital's rule is to be applied twice. 
However this approach is not adopted here. In general, 
the density at the surface is very high, particularly for 
the cold wall case. It has a very larg? gradient at the 
body surface. So density at the surface is calculated from 
the equation of state i.e. equation (4.47) while the pressure 
at the surface is calculated from the expression. 


(4.74) 


where p_ is calculated from equation (4.48) by evaluating it 


at the surface. This gives. 
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(4.75) 
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This process of iteration continues till proper convergence 
of the variables is achieved. The value of n^^ is corrected 
in each iteration till the temperature profile reaches the 
outer boundary conditions asymptotically. It is found that 
at this stage, other variables of the flow also approach the 
freestream values asymptotically. 


4.5 Grid Adaptation 

Ip. the merged layer flow gradients of the flow variables 
are large near the surface and in the shock wave like outer 
region of the flow field whose location is not known a priori. 
Moreover the shock wave like region develops as the integration 
proceeds. Under such flow conditions the adaptive grid technique 
has to be employed so that more grid points move automatically 
to the region of high gradients of the flow field. In the 
present investigation the following weight function is used 


w 


'Ei> 


= 1 + A 


df 


d^f. 


mi +B l^j _ for 1 ^ i ^ NDIV+l 

(4.76) 


which is averaged as 

w( = 


w(Ci) + 


for 1 i ^ NDIV 


Here, f is any dependent variable and A and B are constants 
which indicate the relative importance of the gradients 
and the curvature of the profile of the variable f . Following 
Nakahashi and Deiwert^^, the grid size in the physical plane is 
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related to the weight function in the computational 
plane in the following manner: 


A 


w (Cj_) 


NDIV 

E 


i =1 


1 

w(5^) 


(4.77) 


where n^^ = EAr^. 

It could be seen from equation (4.77) that as the 
weight function increases in magnitude, the grid size 
decreases and vice versa. 

Proper choice of the variable f in equation (4.76) 
is made by considering the conditions at 92 Km. altitude of 
the shuttle flight conditions using the date from Table (4.1) 
Here, Re^ is calculated from Sutherland viscosity law. 

This choice of f is not altered for other altitudes. 

Figure (4.2) could be referred to understand the following 
discussion more easily, 

4.5.1 Choice of the Constants A and B 

There is no fixed criteria as how to choose the 
constants A and B. The best and obvious method, for the 
present, is trial and error. The dependent variable in the 
weight function is chosen either as temperature alone or 
the normal component of velocity alone or a combination 
of both in the manner stated in the next section. For 
each of these choices the value of A is varied from 0.2 
to 1.6 in steps of 0.2. For each value of A, the following 
set of values of B are tried : 0.0, 0.15, 0.20, 0.4 and 


0 . 6 . 




FIG.4.2 PROFILES OF Uo,-Vo,To AND Po WITH UNIFORM 
GRIDS 
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The best choice of A ana B is the one which leads 
to faster convergence, resolves the flow field adequately 
and gives accurate values of at the surface. For the 
prescribed conditions at 92 Kiji, it was found that the values 
A = 1 and B = 0.15 satisfied best the above criterion. 

The gradients and curvatures in the shoch wave 
like region increase as the Reynolds number inreases . A 
single set of values of A and B is not adeguate to resolve 
the flow in the entire range of the Reynolds number. It 
was found that A = 1 and 0.0 £ B < 0.4 gives satisfactory 
results for the entire range of Reynolds number tested. 

It was further found that the constants A .and B also 
depend upon the number of grid points chosen .This point is 
illustrated at appropriate places in the section entitled 
"Discussion of Results". 

4.5.2 Choice of the Dependent Variable in the Weight 

Function w 

At low Reynolds number, the shock is diffused 
and the gradients are large near the surface than in the 
shock wave like outer region (See Fig. 4.3) . At higher 
values of Reynolds number, a distinct shock is found and 
the gradients of the flow variables and v^ are large 
in the shock wave like region in comparison to the gradients 
in the inner viscous ' region near the surface. In 

the shock wave like region, the curvature of the profiles 
of the normal component of velocity is greater than the 
curvature of the profile of the temperature (See Figs. 4. 2 
and 4.7) , Thus, it was found that when temperature alone 


was taken as the dependent variable in the weight function, 
much finer grids were obtained near the surface and 
coarser grids in the outer shock wave like region. When 
the normal component of velocity only was used as the depen- 
dent variable in the weight function, finer grids were obtained 
in the outer shock wave like region and coarse grids were 
obtained in the viscous region near the surface. At this 
stage further trials were made with the profiles of tangential 
component of velocity u^, density pressure p^ and mass 

flux Pq^o dependent variable in the weight function, 

under the prescribed conditions and with the set of constants 
A = 1.0 and B = 0.2. It was found that each of them gave 
unsatisfactory grids. As such, it can safely be concluded 
from the above that temperature ought to be used as the 
dependent variable in the weight function to resolve the flow 
near the surface and the normal component of velocity as 
the dependent variable to resolve the flow field in the rest 
of the field. From a comparative study of the gradients and 
curvatures of the profiles of temperature and normal component 
of velocity in most of the cases, it was found desirable 
to define the weight function with reference to the temperature 
for the first ten grid points and with reference to the normal 
component of velocity for the rest of the grid points in 


the flow field. 
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4.6 Discussion of Results ; 

Two kinds of study were carried out to determine the 
effectiveness and usefulness of the adaptive grid technique. 

The first study was carried out (1) at = 24.95 and Re =27.74 
with cold wall condition and (2), at = 10.0 and Re^ = 87.08 
with insulated wall conditions. 

For cold wall case, using linear profile for u^, and 
p^and quadratic profile for v^ as the initial profiles, 
converged solutions are sought, on both uniform and adapted 
grids, by reducing the number of grid points from 61 to 16. 
Coniputations are carried out with the constant B taking the 
values of 0.0 and 0.15 respectively. When B = 0 curvature term 
in w is not considered and when B = 0.15 curvature term is 
considered. CPU time required and number of iterations taken 
to arrive at the converged results are given in Table (4.2) 
along with the corresponding Cjj values for 61, 31 and 16 grids. 

It could be seen from Table (4.2) that when the 
number of grid points is 61 or 31, solutions with the 
adapted grids took considerably more number of iterations and 
larger CPU-time in comparison to the solution with uni fojrm 
grids. Both the approaches gave almost the same value 
of C„ and the maximum temperature (Fig. 4.4). Results 

; ' II ■ . ■ ■ 

with the adapted grids with B = 0.15 took larger number of 
iterations and more CPU-time. It could be seen from 
Table (4.2) that solutions with uniform grids failed to 
converge wfth 16 grid points and solutions with adapted grid 
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with B = 0 gave a value of C^, which is much different 
from the value of 0.433 obtained with 61 grid points on 
a uniform grid. But with B = 0.15, the resulting value 
of from solutions on the adapted grid agreed well with 
the value obtained on a uniform grid with 61 grid points, 
although compared to the case with B = 0.0 and 16 grid 
points it took more time and iterations. Figure (4.3) 
shows the profiles on a uniform grid with 61 grid points 
and an adapted with 16 grid points and with B = 0.15. 

Both the profiles agree well with each other. 

Table 4.3 gives the results for . the condition of 
= 10..0 and Re^ = 87.08 and using insulated wall boundary 
condition. Using 61 grid points, it could be seen that 
the results with adapted grid with B = 0.0 in the weight 
function took lesser number of iterations and lesser CPU-time 
than the results with uniform grids . These conditions are 
quite different from the results shown in Table 4.2 for 
the shuttle flight conditions at 99.4.9 Km altitude and 
cold wall conditions. The reason for the anamoly lies 
in the vastly different nature of the profiles of temperature 
and normal component of velocity between the two cases. 

For example tiv. W'-, for the insulated wall case, the 
temperature profile in the viscous region near the wall is 
almost perpendicular to the wall and varies almost linearly 
in the outer, thick shock wave like region. Thus the variation 




(a) TEMPERATURE PROFILES (b) TANGENTIAL VELOCITY PROFILES 

FIG. A. 3 COMPARISON OF THE FLOW FIELD WITH UNIFORM AND ADAPTED GRIDS 
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in tlie adapted grid sizes due to variation in curvature 
do not play any significant role. Simply gradient 
variations are enough to give rsasonarly accurate grid 
size for faster convergence. 

In the second study, computations of the governing 
ec'ua cions -.or the merged layer have been carried out for 
the shuttle flight condition from 104.93 Km. to 85.74 
Km altitude (Table 4,1). Here the Reynolds number was 
calculated using Sutherland law of viscosity and it varied 
frcm 10.28 to 260 while the Mach number varied slightly in 
the entire trajectory. Here, n^^ is not allowed to increase 
beyond its value at Re^ = 110. Variation of the 
maximum temperature with uniform and adapted grids 

against Reynolds number is plotted in Fig. 4.4. Also, the 
value of post shock temperature obtained from the Rankine- 
Hugoniot relations is plotted in Fig. 4.4. From Fig. 4.4, 
it could be seen that the values of f is exactly the 

Tn3X 

same for both uniform and adapted grids upto Re^ = 160. 

Also, upto Re^ = 160, the detailed structure of the flow 
for the two approaches agreed reasonably well (see Fig. 4.3 
for Re.= 24.73). For Re_ above 70, a distinct shock and 

O u 

a boundary layer with an invircid flow in between is -Formed. 
As such the predicted f_ „ agreed with R.H. Value of 0.972 

witli reasQiiiable accuracy/ maximum error being 1*5 percent. 

For tbe range of Reynolds number 160 «< Re^ < 260/ the 

results with uniform grids ' did -converge by satisfying the 
convergence criterion but they, -are .in error in that de- 

creasea with increase of Reynolds number which is physically 


max. 



FIG *•* 


iMPARISON OF '"'f” UNIFORM 

,apteo grids with rankine-humniot post-shock 

:mpeRATURE for 10-28S Re^S *30 



unplausible. On the other hand, the results with adapted 

grids for 160 <, Re^ £ 260 gave the correct variation of 

f as it agreed with R.H values. Besides, with adapted 
max 

grids the computation could be extended upto Re^ = 430. 

At this value of Re it could be expected reasonably well that 

o 

the results of the merged layer theory has approached the 
upper limit of the boundary layer theory- Hence, from the 
above discussion it could be concluded that in the denser 
region of the shuttle flight, adapted grid solution gave 
physically accurate results and extended the range of 
■yQ5_id.ity of the merged layer model of the flow field to a 
reasonably large Reynolds number. 


In Fig. 4.5 the converged profiles of u^ ,-v^ and 
on both uniform and adapted grids at Re^ = 230.08 are plotted. 

In Fig. 4.6 the temperature profile alone is plotted for 
both uniform and adapted grids at Re^ = 260.08. It could 
be seen from Fig. 4.6 with uniform grids, the temperature 
tends to decrease in the inviscid zone and attains a value 

n-F T = 0 89 while with adapted grid, the temperature 

Max 

almost remains constant ana attains a maximum value ot 
0.98. Moreover the number of grid points in the shook 
wave like region on uniform grids is thtee whic 
hardly sufficient to resolve such steep gradients as those 
existing in the shock. But the number of grid points in 

the shook wave like region on the adapted grids is 9. Considering 
that the total number of grid points is only 31,the ability of the 




FIG. 4.5 COMPARISON 


OF Oo f Vo AND To AT Re© « 230 




RG. 4.6 temperature PROFILES WITH UNIFORM AND 
ADAPTED GRIDS AT Reo=260 
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adapting technique to put 9 grid points in the shock wave 
like region is highly desirable. In Fig. 4.7, the profiles 
of u^, “^o* ^o ^o adapted grid are shown at 

Re = 400. These profiles indicate a very thin and strong 
shock. 

In Figs. 4.8 and 4.9, heat transfer coefficient at 

the surface C is plotted against Knudsen's number and 
H 

Reynold ' s number respectively . 




FIG. 4 8 COMPARISON OF HEAT TRANSFER COEFICIENT , Ch ,W!TH UNIFORM 
AND ADAPTED GRIDS FROM 104.93KMTO 85-74 KM ALTITUDE 
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Table 4.1 : Input Data Shuttle Flight Conditions 


Altitude 

Km 

Body 

radius 

cm 

'^wall 

°K 

Po, 

00 

gm/cm^ 

T 

00 

U 

00 

cpd/r 

o 

M 

104.93 

136.2 

506.5 

0,246x10“^ 

223.0 

747000 

10.28 

24.95 

99.49 

136.2 

705.0 

0.591x10"^ 

190.0 

750000 

34.74 

27.14 

92.35 

129.6 

921.5 

0.218x10"® 

180.0 

750000 

87.08 

27-88 

85.74 

132.2 

1197.0 

0.637x10"® 

198.9 

753000 

258.08 

26.63 
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Table 4.2 : Shuttle Flight Conditions - Cold Wall 

Alt = 99.49 Km, M = 27.14 Re^ = 24.74, = 1.4, T = 0.025 

00 QOO ' • l/y* 


Points 

Nature of 
grid 

the 

CPU-time 
min: sec 

ITER 

C 

H 

61 

Uniform 


1:57.63 

1600 

0.433 

61 


0 

2:15.18 

1800 

0.444 

61 

^asp{?ag , 

0.15 

3:27.74 

2300 

0.433 

31 

Uniform 


40.57 

1100 

0.401 

31 

Adapted 
A=l» B = 0 


45.58 

1200 

0.426 

31 

Adapted 
A=l; B = 0 

.15 

1:47.75 

2000 

0.414 

16 

Uniform 

did 

not converge 

upto 4000 

iterations 

16 

Adapted 

A = 1 ; B=6 

- 

24.24 

1200 

0 . 365 

16 

Adapted 

A=1;B=0,15 


22.48 

1300 

0.428 



Table 4,3 : Insulated Wall 

= 10,00; Re^^ = 87,08 ;y = 1.667; Pr = 0.75; = 1.0 


Grid points 


Nature of the 
grid 


CPU-time 

min/sec 


ITER 


'H 


61 

61 

21 


Uniform 1:28.84 1400 -0.006 

Adapted 1:15.57 1200 -0.005 

A=1;B=0 

Uniform did not converge upto 4000 iterations 
Adapted 


21 


20.32 


800 


+ 0.002 
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